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1  Introduction 

Swimming  organisms  adjust  their  behavior  in  response  to  environmental  conditions  and  form 
structures  such  as  patches,  swarms,  and  schools.  Disadvantages  of  these  tendencies  include 
tougher  competition  for  food,  the  attraction  of  predators;  some  hydrodynamic  disadvantages 
-  turbulent  wakes  disrupt  movement,  and  more  energy  is  required  to  swim,  etc.  There  are 
several  advantages  however,  such  as  enhanced  reproduction,  predator  avoidance,  and  the 
easier  search  for  prey;  some  hydrodynamic  advantages  -  extra  turbulence  brings  higher 
encounter  rates,  and  coordinated  swimming.  Here,  we  study  a  model  for  such  organisms’ 
behavior. 

2  Lagrangian  Dynamics 

In  order  to  consider  the  model  mentioned  above,  we  start  with  the  Lagrangian  dynamics  of 
an  individual  or  molecule  of  fluid.  The  stochastic  differential  equations  for  the  position  X 
and  velocity  U  are 

dXi  =  Uidt ,  (1) 

dUi  =  Aidt  +  PijdWj,  (2) 

where  A  is  the  acceleration  produced  by  deterministic,  or  large-scale,  forces,  flijdWj  is  the 
random  acceleration,  with  the  random  increment  dWj  satisfying  (dWi)  =  0  and  (dW{  dW ))  — 
Sijdt ,  and  flij  will  be  a  diagonal  matrix.  The  equations  (1)  and  (2)  also  hold  for  continuous 
X  and  discrete  U.  As  an  example,  consider  a  drag  law  for  the  acceleration 

A  =  — r(U  —  u),  (3) 

with  u  being  the  water  velocity  and  /3ij  =  fdSij.  Then  the  dispersion  is  determined  by  /?  and 
r;  from  the  equations,  by  assuming  that  the  statistics  become  time-independent  for  large  t, 
( Ui(t ))  =  (Ui(t  +  dt))  in  the  long-time  limit  or  solving  the  equations  (1)  and  (2)  by  utilizing 
Ito’s  lemma  (see,  for  example,  [2]),  we  can  show  that 


(Ui)  Ui , 

(4) 

(5) 

(32 

(XiWXjit))  ->  (V(O)V(O))  +  ^ Sijt . 

(6) 
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Figure  2:  Histogram  of  velocities  from  numerical  simulations  of  Lagrangian  particles  with 
the  parameters  (a)  r  =  l,/3  =  0.5, u  —  0.5, t  =  30  (b)  r  =  1, /3  =  0.5,?/  0.25, £  300  (c) 

r  =  2,/3  =  0.5,  ?/  =  0.5,  t  =  30.  The  blue  and  the  green  lines  in  each  panel  represent  the 
velocities  in  the  x-  and  y-directions,  respectively. 
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As  diffusion  coefficient  k,  is  k  =  ( dXidXj )  /(2d£),  the  latter  corresponding  to  a  diffusivity 
of  k  —  /32/2r 2,  which  means  that  the  area  grows  like  4 hit  and  the  velocity  variance  is  rn. 

Figure  1  contains  snapshots  of  the  distribution  of  particles  for  different  combinations  of 
the  parameters  r  and  /3.  We  can  see  that  the  area  occupied  by  the  particles  grows  more 
for  larger  (3  (see  Fig.  1(b)  and  (c))  and  smaller  r(see  Fig.l  (b)  and  (d)).  This  can  also 
be  confirmed  by  the  probability  density  functions  (not  normalized  to  one)  of  the  x-  and 
y- velocities  shown  in  Fig. 2. 

3  Grouping  Mechanisms 

To  consider  a  variety  of  motions,  we  use  the  form 

dXi  =  Uidt , 

dU i  =  —r(Ui  -m-  Vx)dt  +  fajdWj.  (7) 

Here,  V  denotes  the  preferred  swimming  velocity,  swimming  towards  the  surface  for  exam¬ 
ple.  Now  we  distinguish  between  various  behaviors  of  swimming  organisms.  The  properties 
described  below  maybe  used  individually  or  together,  depending  upon  the  organisms  con¬ 
sidered. 

•  taxis:  V  depends  on  gradient  of  cue  field  VC(x,t)  (Fig. 3),  where  cue  C(x,  t)  may  be 
environmental  (food,  light,  depth,  etc.)  or  social  (positions  of  neighbors,  etc.).  This 
describes  a  large-scale  preferred  velocity  that  the  group  tends  to. 


Figure  3:  Representation  of  taxis. 

Here,  we  define  taxis  as  a  preference  for  moving  up  the  gradient  of  the  cue  field, 

V  =  aVC(x). 


As  an  example,  we  consider 


C  —  Cq  [1  —  cos(fcx)]  /2  (8) 

so  that  V  =  Vo  sin (kx)x  and  A  =  —  r  [U  —  u  —  Vo  sin(fcx)x].  Figures  4  and  5  show  the 
distribution  and  the  PDF  (not  normalized  to  one)  of  the  organisms  taken  from  numerical 


47 


simulations,  respectively.  Taxis  on  a  spatially  and  temporally  fixed  cue  field  has  its  center 
of  aggregation  where  V  •  Vprey  (or  V2(7)  is  most  negative  when  there  is  no  advection  (see 
Fig. 4  (b)),  and  the  advection  can  shift  this  downstream  (see  Fig. 4  (b)-(d)),  which  decreases 
the  strength  of  the  aggregation  as  a  result  (see  Fig. 4  (b),  (c),  and  (d)).  We  can  also  see 
that  the  diffusion  coefficient  k  —  /32 /2r2  controls  the  width  of  the  aggregation  (Figs. 4,  5). 


t  =  0.000  u=0.  0=0.7,  t=30. 000  u  =  0.5.  0=0.7,  t=30.000 


X 


X 


Figure  4:  Taxis:  Numerical  simulation  results  of  the  example  described  in  Section  3  showing 
organism  distribution  (blue  dots)  and  the  profile  of  C  (green  line)  with  (a)  the  initial 
conditions  for  all  cases,  (b)  u  =  0,/3  =  0.7,  t  =  30,  (c)  u  =  0.5 ,/3  =  0.7,  t  =  30,  (d) 
u  =  0.5,  0  =  0.35,  t  =  30,  (e)  u  =  1.1,0  =  0.7,  t  =  30. 


We  then  define  social  taxis  as  grouping  when  the  cue  field  is  defined  in  terms  of  the 
positions  of  other  organisms; 


c(x)  =  x; 


-IX'  -XI4- -IX'  -XI2 
4  1  21  1 


|X'-X|  <  1, 


(9) 


so  that  the  VC  can  be  expressed  as  VC  =  (X  —  X,)io(|(X  —  X')|)  with  weight  function  w 
having  the  profile  shown  in  Fig. 6.  Using  this  form  of  social  taxis,  the  temporal  variation  of 
the  aggregation  of  the  considered  organisms  depends  upon  the  displacement  of  each  set  of 
two  organisms,  i.e.  |X7  —  X|,  and  so  groups  close  to  each  other  merge  as  time  progresses  as 
shown  in  Fig. 7  (note  the  groups  in  red  circles). 


•  kinesis:  0  depends  on  the  cue  field: 

/?  =  /3(C).  (10) 

This  describes  an  individual’s  tendency  to  move  randomly,  when  they,  for  example, 
find  some  food  close  to  them  (Fig. 8),  and  is  a  more  primitive  response  in  which  the 
random  accelerations  increase  or  decrease  depending  on  the  cue  field. 
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Figure  5:  Taxis:  Histograms  from  numerical  simulations  of  organism  distribution  in  the 
x-direction  for  (a)  u  =  0,/3  =  0.7,  t  =  30,  (b)  u  —  0.5 ,/3  =  0.7,  t  =  30,  (c)  u  —  0.5,  f3  = 
0.35,  t  =  30,  (d)  u  =  !.!,/?  =  0.7, t  =  30. 


Figure  6:  An  example  of  the  weight  function  w( |X'  —  X|)  used  in  social  taxis.  Here, 
w(x)  =  ( \x\ 3  —  |x|)/x. 
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Figure  7:  Social  Taxis:  The  temporal  variation  of  organism  distributions  from  numerical 
simulations.  From  left  to  right,  the  panels  show  the  organism  distribution  at  times  t  — 
0,  t  =  27,  and  t  =  28. 
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Figure  8:  Representation  of  kinesis. 

Let  /3  =  /3o  —  /?iCyCb,  for  example,  and  consider  C  given  by  (8)  again.  Figures  9  and  10 
show  example  distributions  of  organisms  from  numerical  simulations.  Kinesis  can  produce 
aggregation  since  /?  decreases  as  C  increases,  which  means  dU/dt  decreases  as  C  increases. 
The  groups  tend  to  be  looser  however,  depending  upon  flmax/ flmim  and  u  as  mentioned 
before. 

As  another  example  of  kinesis,  we  assume  that  the  organisms  turn  more  frequently  and 
have  smaller  mean  free  path  in  the  presence  of  many  neighbors.  Now  we  consider  /?  in  the 
form  fd  oc  exp(— C/Co),  and  write  social  kinesis  as 

C(X)  =  Yi  [1  -  |x-x'|2]  ,  |X-X'|<1. 

X' 

This  does  not  include  the  fourth  order  term  as  in  (9),  since  social  kinesis  appears  in  (3 
as  a  non-derivative  form.  The  possible  situations  that  can  be  produced  are  some  random 
turning,  more  random  turning,  and  the  avoidance  of  others  as  shown  in  Fig.ll. 
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Figure  9:  Kinesis:  Distributions  of  organisms  from  numerical  simulations  with  (a)  initial 
time  for  all  cases,  (b)  u  =  0,/?o  =  2.7, /?i  =  2.0,  t  —  30,  (c)  u  =  0,/3q  =  2.7,  —  2.6,  t  —  30, 

(d)  u  =  0.5,  =  2.7,  ft  =  2.4,  t  =  30. 


(a)  u  =  0.  0o=  2.1,  fa-2,  t  =  30.  000  (b)  u  =  0.  0o=  Z  7,  0,  =  Z  6.  t  =  30.  000  (c)  u  =  0.  5.  0o=  Z  7,  0,  =  Z  4,  t  =  30.  000 


Figure  10:  Kinesis:  Histograms  of  organism  distributions  in  the  x-direction  from  numerical 
simulations  for  (a)  u  =  0,/?o  =  2.7,  ft  2.0,  £  =  30,  (b)  u  =  0, /So  =  2.7,  /?i  =  2.6,  t  =  30, 
(c)  7/  =  0.5,  /3o  =  2.7,  ft  =  2.4,  t  =  30. 
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Time 


Figure  11:  Social  kinesis,  (a)-(d)  Distributions  of  organisms  from  numerical  simulations 
showing  (a)  the  initial  distribution  for  all  cases,  distributions  after  the  time  shown  in  the 
cases  (b)  with  some  random  turning,  (c)  more  frequent  random  turning,  and  (d)  organisms 
avoiding  neighbors;  and  (e)  the  temporal  variation  of  spatial  average  of  C  for  panels  (c) 
(blue)  and  (d)  (black). 


52 


In  the  case  that  organisms  turn  randomly  occasionally,  the  system  ends  up  with  a 
large  group  of  organisms  (Fig.  11(b)),  whilst  the  other  two  cases  end  up  with  much  looser 
groups  (Fig. 11(c)  and  (d)).  Comparing  the  cases  of  turning  randomly  and  avoidance,  the 
avoidance  case  shows  stronger  grouping,  since  the  small  mean  free  path  interferes  with 
the  avoiding  behaviour.  This  can  also  be  verified  by  the  temporal  variation  of  the  spatial 
average  of  the  social  kinesis  (7,  which  measures  the  strength  of  the  grouping;  the  avoidance 
case  shows  a  much  larger  value  of  average  C  than  that  of  the  case  with  frequent  random 
turning(Fig.  1 1  (e) ) . 

•  Schooling:  V  depends  on  neighbor’s  U  with  |V|  having  a  fixed  value  (Fig. 12):  V  = 
V (Uneighbors)  •  This  describes  the  behavior  of  the  organisms  that  tend  to  swim  simi¬ 
larly  to  their  neighbors. 


■p  /  „  p 
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p  pp  fpp 
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P  P 
PP 


Figure  12:  Representation  of  schooling. 

The  preferred  direction  of  the  swimming  organisms  results  from  a  combination  of  at¬ 
traction  and  alignment  tendencies,  and  so  schooling  can  be  represented  as 

v  =  VbVi/|Vi|, 

Vi  =  a  yyx'  -  X)w(|X'  -  X|)  +  U'w(|X'  -  X|). 

X'  X' 

The  first  term  represents  attraction  to  neighbors,  a  controls  the  strength  of  the  attraction; 
larger  a  gives  stronger  schooling  as  seen  in  Figure  13  which  shows  the  temporal  variation  of 
organism  distributions  from  numerical  simulations  for  different  a.  Paying  attention  to  the 
groups  surrounded  by  circles,  we  notice  that  the  case  a  =  0.3  shows  weak  schooling  and  the 
groups  tend  to  break  apart  easily,  whilst  the  case  a  —  0.7  shows  much  stronger  schooling. 
The  second  term  represents  the  tendency  of  organisms  to  align  their  swimming  with  their 
neighbors.  The  distance  between  two  organisms  at  which  alignment  becomes  effective  tends 
to  be  smaller  than  that  for  attraction,  as  shown  in  Fig.  14.  Here,  the  choice  of  the  weighting 
function  w  affects  the  internal  structure  of  self-organized  groups. 
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Figure  13:  Schooling:  Snapshots  of  organism  distributions  from  numerical  simulations,  (a) 
the  initial  condition  for  all  cases,  (b-1)  a  —  0.7,  t  =  580,  (b-2)  a  —  0.7,  t  —  585,  (b-3) 
a  =  0.7,  t  =  590,  (c-1)  a  =  0.3,  t  =  615,  (c-2)  a  =  0.3,  t  =  620,  (c-3)  a  =  0.3,  t  =  625,  (c-4) 
a  0.3,  t  630. 


Figure  14:  Schooling:  Areas  of  attraction  (outer  circle)  and  alignment  (inner  circle).  The 
gray  dots  represent  organisms. 
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The  Boltzmann  Equation 


We  are  interested  in  finding  the  form  of  the  probability  density  function  of  the  system, 
P(Xi,  Ui,  X2,  U2, ,  t ).  The  evolution  of  this  quantity  is  described  by 

dv  r 

—  =  - VXi  •  (UiV)  -  VUj  •  (A -P)  +  J  dnU' T(U  -  U'\x,  U',  t )  V.  (11) 

If  V  is  sharply  peaked  in  velocity,  the  last  term  on  the  RHS,  which  can  be  thought  of 
as  a  ‘collision’  rate,  can  alternatively  be  written  as 

<92  f 

SUiVj  V  2  )' 

If  we  assume  independence  of  organisms  i,  which  is  not  true  for  social  behavior,  the 
probabilities  are  equal  for  all  i,  and  we  can  rewrite  the  evolution  as 

^  =  -Vx  •  (UP)  -  Vu  •  (AV)  +  \vl(32v-  (12) 

The  first  term  on  the  RHS  represents  advection  in  the  x-direction,  proportional  to 
the  velocity  at  that  point,  the  second  term  represents  advection  towards  the  x-axis  (if 
A  ^  —  rU),  and  the  third  term  represents  diffusion.  The  first  two  terms  are  depicted  in 
Figure  15.  As  can  be  seen,  organisms  will  be  directed  towards  the  x-axis  and  parallel  to 
it,  resulting  formation  of  a  thin  filament,  which  is  restrained  from  reaching  zero  width  by 
diffusion.  In  taxis ,  A,  the  advection  towards  the  x-axis,  depends  on  position  X,  and  in 
kinesis ,  /?,  the  diffusive  term,  depends  on  position. 


Figure  15:  Representation  of  the  directions  of  the  first  (dashed)  and  second  (dotted)  terms 
on  the  RHS  of  equation  (12).  In  effect,  organisms  filament  along  the  x-axis,  although  the 
diffusive  term  stops  the  width  from  going  to  zero. 
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Examples 


Example  solutions  to  equation  (12),  if  we  define  as  before  /3  =  /3q  —  PiC/Cq  and  A  = 
— r(U  —  u),  can  be  seen  in  Figure  16.  Figure  16b)  shows  that,  with  diffusion  constant  in 
space,  i.e.  /? i  =  0,  the  organisms  group  in  velocity,  but  not  in  space.  Adding  a  background 
flow,  u,  or  diffusion  dependent  on  the  cue  field  C(x),  causes  grouping  in  both  |U|  and  |X|, 
see  Figures  16d)  and  g). 


Figure  16:  Plots  showing  the  results  of  solving  the  Boltzmann  equation  (12).  In  a)  c)  and  f), 
the  velocity  of  the  organisms,  U,  is  plotted  vs.  |U|  and  |X|,  in  b)  d)  and  g),  the  probability 
density  function  V  at  large  t  is  plotted  vs.  |U|  and  |X|  and  in  e)  and  h),  the  density  of 
organisms  p  at  large  t  is  plotted  vs.  |X|.  Plots  a)  and  b)  have  /3q  =  2,  (3\  =  0  and  u  =  0. 
Plots  c)  d)  and  e)  have  /3q  =  1,  =  0  and  u  =  2.  Plots  f)  g)  and  h)  have  /3o  =  2.7,  /A  =  2.4 

and  u  =  0.  All  quantities  have  been  non-dimensionalised. 


Density 

Let  us  define  the  density,  p  —  f  dXJV,  and  thus  the  evolution  is  defined 

^  vx  J  dXJ  UV  =  -VxUft  (13) 
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where  the  bar  indicates  a  space-mean. 


Velocity 

This  leads  to  the  following  expression  for  the  evolution  of  the  quantity  pU^: 

=  -VXiOTu7  +  p\Wj)-  (14) 

We  can  also  write  —VU(AV)  =  —  Vu/(— rU  —  rXJ')V ,  and  then  replace  these  terms  back 
into  equation  (12).  It  can  be  seen  that  the  system  is  not  closed,  but  if  we  assume  that  the 
three  terms  in  (12)  roughly  balance,  i.e. 

_ _ B2 

U^Uj  «  XJj  + 

then  the  system  is  closed,  and  V  is  now  quasi-Gaussian.  We  can  thus  derive  the  closed  mass 
equation 

^  +  V-(Up)  =  0,  (15) 

and  the  momentum  equation 

dtJ-  _  _  IB2  _ 

+  Uj-VjUi  +  =  — r(U;  -  iii  -  V).  (16) 

Drag  Dominated  Case 

If  drag  dominates,  then  we  can  set  the  velocity  U,  and  so  the  momentum  equation  becomes 

dp  1  B 2 

-£+V-(uP  +  VP)  =  V--V^P,  (17) 

Taxis  vs.  Kinesis 

Let  us  re-write  equation  (17)  in  terms  of  diffusivities: 

^  +  V  •  (u  +  V  —  -V(rzt)p)  =  V  •  /«Vp,  (18) 

at  r 

where  the  diffusivity  k  =  /32/2r2  and  is  spatially  varying.  The  third  term  on  the  LHS 
represents  taxis ,  and  the  fourth  kinesis .  In  this  form  it  can  be  seen  that  both  act  similarly 
on  the  system,  converging  velocities  to  certain  regions,  either  those  with  smaller  velocities 
(taxis)  or  low  diffusivities  (kinesis).  This  can  also  be  seen  in  Figure  17,  which  shows  the 
initial  and  final  states  of  organisms  experiencing  either  (taxis)  or  (kinesis),  as  well  as  the 
final  density  states  in  both  cases.  However,  kinesis  has  similar  diffusivities  and  velocities, 
i.e.  Peclet  number  is  of  order  1,  while  the  velocity  in  taxis  can  be  higher. 
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Figure  17:  Numerical  simulations  of  the  system,  assuming  either  taxis  in  a)-c),  i.e.  V  = 
aVC  is  spatially  varying  but  /?  is  not,  or  kinesis  in  d)-f),  i.e.  /3  is  spatially  varying  but  V  is 
not.  The  green  line  in  Panels  a)  and  b)  shows  the  distribution  of  the  concentration  field  C. 
Panels  a),  d)  and  b),  e)  show  the  initial  and  final  distributions  of  organisms.  Panels  c)  and 
f)  show  the  solutions  of  the  Boltzmann  equation  (12),  [blue  line],  the  density  equation  (13) 
[green  line]  and  a  histogram  of  the  density  in  the  final  snapshot  of  the  simulation.  Note  that 
the  blue  and  green  lines  lie  on  top  of  each  other  in  Panel  f).  These  show  that,  while  the 
biological  mechanisms  behind  them  is  different,  both  processes  produce  similar  final  density 
distributions. 
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Stability  with  social  behavior 

In  the  absence  of  a  background  flow,  i.e.  u  =  0,  we  redefine  V  =  V  +  ^  V  •  (Yk),  then 

^  +  v  •  (Vp)  =  V  •  KVp.  (19) 

The  social  behavior  is  specified  by  making  k,  a  functional  of  the  density,  defined  in  terms  of 
a  weighting  function  w,  such  that 


k  =  K( 


ds  w(s)p(K  +  s)), 


and  w  takes  the  form  of  decreasing  with  s,  which  represents  the  sensing  distance  of  the 
organism,  see  Figure  6.  If  the  basic  state  density  is  uniform,  then  k,  is  constant  in  space, 
and  we  can  rewrite  p  and  k  as  a  sum  of  the  means  (over-bars)  and  small  perturbations 
(primes),  and  so  the  density  equation  becomes 

f)  Qf 

-Ji  +  pV-V^V-TcVp',  (20) 

where 

V'  =  -Vk'  =  —K'(p)V  J dsw(s)p\X  +  s).  (21) 

Taking  a  Fourier  transform  of  this  equation,  we  can  find  an  expression  for  the  growth 
in  density  perturbations: 

Jlp'  =  [; pK'(p)k2w(k )  -  Kk2]p.  (22) 

If  the  bracketed  term  is  greater  than  zero,  the  system  experiences  positive  growth  see 
Figure  18,  which  shows  the  relationship  between  growth  rate  and  wave  number  k.  In  fact, 
the  system  experiences  explosive  (faster  than  exponential)  growth,  with  the  amplitude  going 
as 

r)A 

—  =  a(k)A+J\fA 3, 

where  N  >  0. 


Evolution  of  Social  Behavior 


I:  Inheritance  of  Traits 


As  before,  we  define  the  velocity  V  as  the  gradient  of  a  scalar,  i.e.  V  =  V0,  see  [7],  where 


=  wq  J dsw(s)Z(X.  +  s),  Jdsw(s)  =  l. 


(23) 


In  this  case,  the  weighting  function  w( s)  takes  the  form  depicted  in  Figure  19  with  repulsion 
at  short  distances. 

We  define  three  different  genetic  types,  or  alleles,  of  the  heterotroph  Z  (which  preys  on 
P),  each  with  two  genes  with  two  possible  genetic  traits,  0  or  1:  Zqq,  Zqi,  Z\\.  Zqo  and  Z\\ 
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Figure  18:  Figure  from  [7],  showing  “the  linear  growth  rate  of  small  perturbations  around 
the  equilibrium  biomass  versus  the  wave  number  of  the  perturbations,  ...  for  different  val¬ 
ues  of  r”,  oc  1/ft,  a  diffusion  time  scale.  “In  a  uniform  distribution  of  particles,  random 
variations  will  spontaneously  lead  to  the  formation  of  patches,  if  there  is  positive  growth 
somewhere  along  the  wave  number  spectrum.  The  gray  line  shows  the  case  where  particle 
motion  is  decoupled  from  biological  activities.  . . .  When  the  diffusion  rate  is  fast  compared 
to  growth  (r  <<  1),  population  dynamics  have  a  negligible  impact  on  the  linear  growth 
rate  of  perturbations  in  the  region  of  instability  (positive  growth).  As  biological  processes 
become  faster,  their  effect  ...  [is  to]  stabilize  the  equilibrium  solution;  the  growth  rate  is 
negative  for  all  wave  numbers  when  r  ~  1” . 


Figure  19:  Figure  from  [7],  showing  the  weighting  function,  w(x),  of  organisms  as  function 
of  their  distance  from  neighbors  located  at  x  =  0,  as  used  in  equation  (23).  uw(x)  = 
exp(—x2/2)  —  exp(—2x2).  Distances  are  expressed  in  units  of  L,  the  characteristic  perception 
length  scale.  The  sensing  radius  has  an  approximate  value  of  2 L.  . . .  Positive  weighting 
indicates  attraction  to  neighbors.  Repulsion  at  short  distance  prevents  over-crowding.” 
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are  recessive,  and  only  Z\\  exhibits  grouping  behavior,  i.e.  Voo  =  Voi  =  0,  Vn  ^  0  .  The 
population  of  each  type  evolves  as  follows: 


—  a9  p  c  —  ~  V  •  (Va^  —  PS7)Zap  +  Map,  (24) 

where  0  <  a  <  1  is  the  conversion  efficiency  of  P  into  biomass  byZ,0<g<lisa  grazing 
rate,  d  is  the  death  rate,  P  is  the  population  of  prey  and  fap  is  the  inheritance  fraction  of 
the  population  Zap.  Ph  is  the  half-saturation  density,  which  provides  a  limitation  to  the 
amount  of  prey  that  can  be  eaten  at  large  P,  i.e.  Z  depends  linearly  on  P  until  the  capacity 
of  Z  to  eat  P  becomes  reached.  The  c  term  represent  the  “Allee”  effect  [1]:  as  the  density  Z 
is  small,  encounters  between  male  and  female  organisms  are  rare  and  depend  quadratically 
on  density,  whereas  for  large  Z  virtually  all  females  have  mated  and  reproduction  depends 
linearly  on  the  population  gets  large,  virtually  all  females  have  mated.  From  a  Mendelian 
table  of  inheritance [4]  we  can  define  the  fractions  fap  as: 


/oo 

/oi 

hi 
Z 2 


Z020  +  Z{){)Z{)i  +  -Z{ 


01  5 


^oo^oi  +  2^00^11  +  x^oi  +  Wn, 


^11  +  Woi  +  ^^01> 

/oo  +  /oi  +  hi  =  (^00  +  Zqi  +  Zii)2 , 


where  this  implies  that,  for  example,  a  new  generation  of  Zqo  is  produced  by  all  interactions 
between  two  Zqo,  half  of  the  interactions  between  Zqo  and  Zqi,  and  1  in  4  of  interactions 
between  two  Z$\. 

We  also  include  a  mutation  term,  Map,  which  allows  for  the  possibility  of  a  population 
reforming  if  it  dies  out  completely: 

Moo  =  -^Zm  +  -p^oi, 

Mqo  —  —  /^oi  +  9Z00  +  iaZu, 

Mn  =  —9Z11  +  -jiZoi, 


where  fi  is  the  fractional  rate  of  mutation.  The  prey,  P  evolve  as  follows; 


dP 

~di 


9 


P  +  PhZ  + 


+  V^VP, 


(25) 


where  b  is  the  growth  rate,  Pc  is  the  half-saturation  density,  and  all  other  quantities  are  as 
before. 


Examples 

There  are  specific  conditions  under  which  the  grouping  type,  11,  can  win  out  with  respect 
to  the  other  two  types.  If  the  grouping  strength  is  not  too  high  and/or  the  diffusivity  is  not 
too  low,  then  it  can  cause  the  extinction  of  the  other  two  species,  due  to  the  reproduction 
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advantages  of  grouping,  see  Figures  20a)  and  b).  However,  if  the  diffusivity  is  too  low  to 
bring  the  grouping  type  sufficient  prey,  then  the  other  two  types  will  dominate,  while  the 
grouping  type  will  die  out  until  the  groups  are  sufficiently  small  to  survive,  see  Figures  20c) 
and  d).  If  the  grouping  species,  11,  is  introduced  at  low  values  to  the  other  two  types,  then 
it  will  slowly  take  over  until  it  reaches  a  critical  fraction,  /c,  at  which  the  other  two  species 
become  extinct  and  the  grouping  type  wins,  see  Figure  21. 


Figure  20:  Final  state  population  distributions  Z^p  from  numerical  simulations  (a)  and  c)) 
shown  next  to  the  fractional  populations  fap  over  time  (b)  and  d)).  In  a)  and  b),  the  species 
Z\\  causes  the  other  two  to  die  out,  due  to  the  advantages  of  grouping  for  reproduction. 
However,  if  the  diffusivity  is  low,  as  in  c)  and  d),  then  the  disadvantages  of  grouping,  i.e. 
the  fast  depletion  of  prey,  mean  that  Zqo  and  Zqi  do  much  better  than  Zn,  which  doesn’t 
completely  die  out. 

If  we  subject  the  system  to  background  stirring  motion,  such  as  Ray  Pierrehumbert’s 
exact  sin/sin  stirring[5],  see  Figure  22(i),  then  we  can  test  whether  grouping  is  still  advan¬ 
tageous  or  not,  see  Figures  22 It  can  be  seen  that,  at  low  grouping  parameters  and 
low  stirring,  the  grouping  type  dies  out,  but  at  large  grouping  parameters  and  large  stirring, 
the  grouping  type  wins.  In  other  regimes,  both  types  survive. 

Thus,  it  is  not  only  the  biological  but  the  dynamical  features  which  can  determine  how 
a  particular  species  will  fare. 
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neutral  evolution - >-  *4 - stable  state- 


Figure  21:  Frequency  of  Zaf 3  versus  time,  showing  the  invasion  of  a  non-grouping  population, 
Zoo,  by  the  grouping  type,  Z\\%  taken  from  [7].  “Time  is  normalized  by  the  mutation 
timescale  . . .  The  new  allele,  introduced  at  low  frequency  at  t  =  0  does  not  affect  fitness 
until  the  frequency  of  grouping  organisms  (type  11)  reaches  /c,  at  which  part  it  starts  to 
aggregate.  ...  In  this  case,  social  behavior  is  a  successful  strategy:  the  grouping  type  wins 
the  competitions  and  the  other  types  become  extinct.  The  ecosystem  is  then  stable  with 
respect  to  mutations  in  the  two-allele  model.” 


II:  Grouping  Parameter 

Examples  of  how  diffusion  and  grouping  affect  the  dynamics  when  limit  cycles  may  be 
present  (through  setting  the  value  of  c  in  (24)),  can  be  seen  in  Figure  23. 

To  find  optimal  parameters,  we  let  populations  with  different  grouping  strengths  com¬ 
pete,  and  note  whether  their  populations  grow  or  decay.  The  effect  of  the  half  saturation 
ratios,  P \  and  c  in  equation  (24),  on  the  optimal  grouping  strength  can  be  seen  in  Figure  24. 

Ill:  Schooling 

As  before,  we  represent  schooling  as  velocity  alignment  with  neighbors  and  attraction  (ex¬ 
cept  very  nearby).  We  add  a  food  supply,  confined  to  the  mixed  layer  at  the  top  of  the 
ocean,  with  variability  produced  by  upwelling  which  brings  further  nutrients.  As  the  mixed 
layer  is  by  definition  assumed  well-mixed,  downwelling  does  not  change  the  concentration 
of  nutrients.  We  specify  a  streamfunction,  and  a  vertical  velocity,  re,  that  depend  on  the 
depth  of  the  thermocline,  h: 

=  jh ,  w  =  ,  (26) 

where  g'  is  the  modified  gravitational  acceleration  and  /  is  the  planetary  vorticity.  We  can 
then  derive  the  evolution  of  the  nutrient,  N : 

dN_  =  w(w>  0)  _  N)  _  ^ NP  +  v.KfluidVN,  (27) 

where  H  is  the  average  depth  of  the  thermocline,  N^eep  is  the  nutrient  value  at  depth,  and 
P  is  the  population  of  the  autotroph,  e.g.  phytoplankton.  Example  biomass  distributions 
from  numerical  simulations  can  be  seen  in  Figure  25. 
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Figure  22:  (i)  depicts  the  form  of  the  Pierrehumberts  exact  sin/sin  stirring,  (ii)  —  ( iv ) 
show  snapshots  from  numerical  simulations  of  Z\\  (left  panels)  and  a  passive  tracer  (right 
panels);  (ii)  and  (Hi)  have  no  diffusion,  i.e.  k  =  0,  (ii)  and  (iv)  are  subject  to  a  medium 
rate  of  stirring,  and  (Hi)  is  not  subject  to  any  stirring,  (iv)  shows  that,  while  diffusion 
eventually  completely  removes  any  gradients  in  the  tracer,  Z±±  still  remains  grouped,  (v) 
shows  a  summary  of  the  conditions  where  grouping  and  non-grouping  Z  survive,  from  [7]. 
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Figure  23:  Numerical  simulations  of  predator  biomass,  Z,  distribution  in  time  and  space, 
from  [6].  “Thin  black  contours  show  the  instantaneous  spatial  distribution,  at  fixed  inter¬ 
vals  in  time.  Thick  black  curve  shows  the  time  series  of  biomass  at  a  fixed  location,  a) 
Without  aggregation  behavior”  (wo  =  0),  “the  predator-prey  system  reaches  a  limit  cy¬ 
cle  after  a  period  of  transients;  the  simulation  shown  here  is  initialized  near  the  unstable 
coexistence  fixed  point,  b)  With  social  behavior”  (wo  —  2),  “and  fast  diffusion,  patches 
form  spontaneously  and  the  system  loses  its  oscillatory  behavior,  c)  With  social  behavior” 
(wq  —  15),  “and  slow  diffusion,  the  system  supports  a  regular  wave  traveling  through  the 
domain;  transients  are  not  shown,  d)  With  social  behavior”  (wo  =  5),  “and  slow  diffusion, 
the  density  field  can  appear  chaotic.” 
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a) 


b) 


Figure  24:  Optimal  strategy  vs  P \  and  c,  from  [6].  x  is  non-dimensionalised  P^,  A  is  non- 
dimensionalised  c,  and  S  is  our  wq.  “Black  symbols  indicate  optimal  parameter  values  of  [wq] 
from  numerical  simulations  of  the  competition  model,  solid  line  is  the  spline  interpolation. 
Gray  symbols  indicate  that  social  behavior  is  not  a  successful  strategy,  in  which  case  non¬ 
grouping  types  . . .  are  the  best  competitors.” 


We  model  foraging  by  assuming  that  the  turning  rate  of  an  organism  is  dependent  on 
the  difference  between  past  feeding,  £  compared  to  current  food  values.  Here  we  consider 
the  turning  rate  59  of  a  heterotroph,  Z,  in  terms  of  the  population  of  its  prey,  P.  The 
memory  of  food  evolves  as  follows: 

=  P(*i)  _  &  ^ 

dt  P(Xi)  +  Ph  r’  1  ’ 

where  r  is  the  memory  period  of  the  heterotroph.  Then  we  model  the  turning  as 

56  =  1  +  tanh (a(&  -  P(X,)),  (29) 

where  a  is  a  constant.  Figures  26  and  27  show  how  the  success  of  the  foraging  Z  depends 
on  how  fast  they  swim.  Foragers  that  are  too  speedy  over-shoot  the  largest  concentrations 
of  prey;  slow  foragers,  however,  don’t  find  the  prey  for  a  long  time.  Figure  28  shows  how 
successful  schooling  is  for  foragers  in  various  levels  of  turbulence  and  resources.  It  can  be 
seen  that,  for  low  levels  of  resources  there  is  no  clear  advantage  to  schooling,  and  it  can 
actually  be  a  disadvantage  at  low  levels  of  cp  the  ratio  of  alignment  to  attraction,  as  defined 
in  equation  (29).  However,  in  cases  with  medium  or  high  resources,  there  are  values  of  a 
that  provide  a  clear  advantage. 

Aggregation 

Observations  of  Antarctic  krill  swarms  can  be  seen  in  Figure  29.  Generally,  during  the  day 
there  tends  to  be  many  small  swarms  near  the  surface,  and  at  night  there  tend  to  be  a 
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a)  Nutrients  (N) 


&}  Phytoplankton  (P) 


Figure  25:  Figure  from  [6],  showing  results  from  numerical  simulations  with  schooling  and 
turbulent  advection.  a)  biomass  of  nutrients  N,  b)  biomass  of  phytoplankton  P .  “Lighter 
shades  of  gray  indicate  larger  values.  The  two  fields  tend  to  be  anti-correlated.” 
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b) 


Figure  26:  Figures  showing  numerical  simulations  of  foraging  Z  with  a)  fast  swimming 
speeds  and  b)  slow  speeds  at  various  non-dimensionalised  times.  The  white  dots  represent 
the  individual  Z,  the  black  tails  represent  the  organism’s  grazing  memory  and  the  color 
scale  represents  the  concentration  of  prey  P.  It  can  be  seen  that  the  fast  swimmers  over¬ 
shoot  and  never  all  end  up  concentrated  on  the  prey,  whereas  the  slow  swimmers  quickly 
group  over  the  high  concentration  of  prey. 
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Food  Intake 


Swi raining  speed 


Figure  27:  Figure  summarizing  the  (non-dimensionalised)  food  intake  of  foraging  Z  with  a 
specific  swimming  speed  from  numerical  simulations.  Clearly  there  is  an  ideal  swimming 
speed  at  which  the  swimmers  find  the  food  within  a  reasonable  time  but  do  not  overshoot 
the  highest  concentration  by  swimming  too  fast. 


few  large  swarms  at  depths  (^250m).  Different  types  of  aggregation  can  be  classified  using 
a  measure  of  polarization:  the  degree  to  which  the  grouping  is  swarm-like  or  school-like. 
Figure  30  shows  results  from  numerical  simulations,  with  the  degree  of  polarization  plotted 
against  the  degree  of  alignment,  a,  and  the  swimming  velocity,  V.  Five  specific  grouping 
phases  are  identified  and  examples  shown  of  each. 


Conclusion 

It  can  be  seen  that  a  wide  range  of  grouping  structures  in  organisms  can  be  simply  mod¬ 
elled.  We  have  seen  that  in  turbulent  environments  it  is  sometimes  advantageous  to  school 
or  group,  depending  on  the  properties  of  the  organisms  themselves  as  well  as  how  the  envi¬ 
ronment  affects  their  supply  of  food.  In  general,  it  seems  that  there  is  the  need  for  a  certain 
amount  of  resource  availability  to  make  it  advantageous  for  foraging  predators  to  school. 
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more  resources 


Figure  28:  The  foraging  success  of  schooling  Z,  taken  from  [6],  “normalized  by  the  success  of 
random  walkers  under  the  same  conditions.  In  each  panel,  the  relative  success  ...  is  shown 
for  different  strategies  (cq  the  ratio  of  alignment  to  attraction”,  see  equation  (29)  “varying 
between  0  and  1).  Each  column  corresponds  to  a  different  level  of  resource  availability;  each 
row  has  different  flow  characteristics.”  In  the  low  resource  panels,  there  is  no  advantage  to 
schooling,  and  in  fact  it  can  be  a  disadvantage  at  low  values  of  a.  However,  with  medium 
to  high  resources,  there  are  values  of  a  that  provide  a  clear  advantage. 
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Figure  29:  Observations  of  Antarctic  krill  from  acoustic  measurements  [3].  Points  represent 
aggregations,  which  are  plotted  as  a  function  of  the  time  of  day  of  observation  (x-axis)  and 
mean  depth  (y-axis).  The  size  and  color  represent  the  mean  density.  The  general  trend 
shows  a  few  patches  of  large  density  at  depth  at  night  and  many  patches  of  small  density 
near  the  surface  during  the  day. 
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Figure  30:  Simulations  results  from  [6].  The  gray  shading  indicates  the  degree  of  po¬ 
larization;  high  vales  (light  gray)  represents  schooled  groups,  low  values  (black)  represent 
swarming  groups.  This  is  plotted  against  cp  degree  of  alignment,  and  V ,  swimming  velocity. 
Labels  I-V  indicate  distinct  phases:  I  milling;  II  swarms;  III  large  schools;  IV  small  schools; 
V  transient  alignment. 
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Lecture  6:  Mixing  in  the  presence  of  sources  and  sinks 


Jean-Luc  Thiffeault 

Notes  by  Sam  Pegler  and  Woosok  Moon 
28  June  2010 


1  Norms 


In  this  section  we  define  a  measure  of  mixing  that  does  not  necessarily  require  diffusion  to 
measure  the  amount  of  homogenization  that  occurs  during  the  mixing  process.  Recall  the 
advection-diffusion  equation 

—  +  u  ■  ve  =  n\72e,  (i) 

Ot 

where  6  is  a  concentration  field  in  a  finite  domain  fi,  with  no-net-flux  boundary  conditions. 
We  assume  without  loss  of  generality  that 


/ 

Jn 


e<m  =  o, 


(2) 


and  define  the  L2-norm.  or  variance,  as 


||0|||=  /  (3) 

Jn 

Recall  from  Lecture  1  that  the  variance  evolves  according  to 

±\\e\\l  =  -2K\\ve\\l  (4) 

and  decays  in  time  as  the  system  mixes.  The  variance  indicates  the  extent  to  which  the 
concentration  has  homogenized  and  is  thus  a  good  measure  of  the  amount  of  mixing  that 
has  occurred.  However,  the  variance  requires  knowledge  of  small  scales  in  0,  which  we  are 
not  necessarily  interested  in.  A  measure  of  how  well-mixed  the  concentration  is  does  not 
necessarily  require  knowledge  of  how  much  homogenization  has  occurred  due  to  diffusion  at 
small  scales.  This  is  more  in  keeping  with  the  definition  of  mixing  in  the  sense  of  ergodic 
theory  [2].  In  this  regard,  we  proceed  to  consider  the  pure  advection  equation 

—  +  u-ve  =  o.  (5) 

Note  that  in  this  case  equation  (4)  predicts  that  the  variance  satisfies 


d 

dt 


(6) 
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